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1.0 BACKGROUND 


This project is part of a larger, ongoing research effort at Cornell University to develop 
theory and simulation methods suitable for predicting crack turning and arrest in integral 
aircraft structures. In addition to the documents included with this report, documentation 
of work completed throughout the entire project, including work accomplished after the 
period of performance for the present program, will be compiled in the form of a doctoral 
dissertation planned for completion by August 2000. 


2.0 ACCOMPLISHMENTS 

A key accomplishment of the present program was the development and implementation of 
the second-order linear-elastic crack turning theory with fracture toughness orthotropy into 
the FRANC2D fracture simulation code. An early version of the code has already found 
use in analyses performed under the Integral Airframe Structure (IAS) Program (Contract 
NASI- 20014, Task 34, ending November, 1998). Correlation of crack path predictions 
with observed trajectories in test specimens for a variety of test conditions was much better 
than would be expected of any other known theory. Since that time, the code has been 
further improved to reflect new T-stress calculation techniques described in Attachment #1 , 
among other improvements. 

In a collaborative effort, the code driving the second-order orthotropic crack turning theory 
was ported to the shell version of FRANC3D to study crack turning and flapping of narrow 
body fuselage structure under the aging aircraft program. 

A theoretical estimate of the characteristic length, r c , associated with crack turning, has been 
developed based on the concept of plastic instability. This parameter is necessary for 
application of the 2 nd order theory. This work is also described in Attachment #1 . 

With regard to the turning criterion, it has been found that a very similar angle prediction to 
the 2nd order theory max stress theory can be obtained based on elastic CTOD theory. The 
solution is sensitive to the T-stress (in the same manner as the maximum tangential stress 
theory), and to the length of the crack extension (which is roughly equivalent to r c ). A 
description of this finding and a proposed extension to the elastic plastic case is described 
in Attachment #2. 

Test plans and specimens were also prepared under the present program for testing under a 
later segment of the overall effort, and an interpolating function for 3D fracture toughness 
orthotropy was also developed to support future development of a 3D crack turning 
simulation capability (applicable to thick structures). 


1 



3.0 ATTACHMENTS 



Evaluation of the Characteristic 
Length and T- stress for 2 nd Order 
Crack Turning Problems 


Richard G. Pettit 
Anthony R. Ingraffea 


Cornell Fracture Group 


Attachment #1 to Final Report 
NASA Contract NAG- 1-20 13 


December, 1999 



ABSTRACT 


In the course of several years of research efforts to predict crack turning and flapping in 
aircraft fuselage structures and other problems related to crack turning, the 2 nd order 
maximum tangential stress theory has been identified as the theory most capable of 
predicting the observed test results. This theory requires knowledge of a material specific 
characteristic length, and also a computation of the stress intensity factors and the T-stress, 
or second order term in the asymptotic stress field in the vicinity of the crack tip. 

A characteristic length, r c , is proposed for ductile materials pertaining to the onset of plastic 
instability, as opposed to the void spacing theories espoused by previous investigators. 

For the plane stress case, an approximate estimate of r c is obtained from the asymptotic 
field for strain hardening materials given by Hutchinson, Rice and Rosengren (HRR). 

A previous study using of high order finite element methods to calculate T-stresses by 
contour integrals resulted in extremely high accuracy values obtained for selected test 
specimen geometries, and a theoretical error estimation parameter was defined. In the 
present study, it is shown that a large portion of the error in finite element computations of 
both K and T are systematic, and can be corrected after the initial solution if the finite 
element implementation utilizes a similar crack tip discretization scheme for all problems. 
This scheme is applied for two-dimensional problems to a both a p-version finite element 
code, showing that sufficiently accurate values of both K, and T can be obtained with fairly 
low order elements if correction is used. T-stress correction coefficients are also developed 
for the singular crack tip rosette utilized in the adaptive mesh finite element code 
FRANC2D, and shown to reduce the error in the computed T-stress significantly. Stress 
intensity factor correction was not attempted for FRANC2D because it employs a highly 
accurate quarter-point scheme to obtain stress intensity factors. 
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1.0 INTRODUCTION 


Crack turning has been recognized as a potentially important crack arrest mechanism for 
pressurized aircraft fuselage structure, and for longitudinal cracks can result in turning and 
flapping as shown in Figure 1 as reported by Maclin [1]. This behavior contains the 
damage, vents the pressure in a controlled manner, and results in obvious damage which 
can be subsequently repaired. Flapping was observed to occur reliably enough during tests 
of thin skinned, relatively narrow-body fuselages that it was utilized as a fail-safe criterion 
on the 707, 727, and 737 fuselages for regions excluding the joint areas. Similar 
phenomena have been observed in unstiffened cylinders by Swift [2], who also reported 
turning and cracking in an experimental fuselage with adhesively bonded stiffeners [3], and 
Pettit [4], who observed crack turning and arrest in integrally stiffened fuselage structure 
with transverse cracks. 

These behaviors were observed in tests, but were never successfully modeled until 
recently. Turning and flapping was observed in narrow-body aircraft, but not wide-body 
aircraft like the 747 aircraft, and no one really understood why. Also, it was found that 
aging aircraft develop multi-site damage, which can potentially alter the crack turning and 
flapping performance [5]. The need for an accurate crack trajectory modeling capability 
was evident. 

In the last decade, a sequence of authors have studied the fuselage flapping phenomenon, 
including Kosai et al [6], Miller et al [7], Potyondy et al [8], Knops [9], and Chen [10], 
Beginning with Potyondy, an adaptive mesh finite element routine was used similar to that 
proposed by Wawrzynek and Ingraffea [11] but extended to three dimensional shell 
problems, which allows the trajectory of the crack to develop naturally in accordance with a 
user-selected crack turning theory. Potyondy used the first order maximum tangential 
stress theory of Erdogan and Sih [12] to predict the crack trajectory of an adhesively 
bonded narrow-body fuselage panel tested at Boeing, as shown in Figure 2. He was able 
to approximate the actual behavior fairly well in the gently curving region until the crack 
grew near to the tear strap, but was unable to predict the path as the crack turned and grew 
parallel to the tear strap, resulting in flapping. 

The work of Kosai, Knops and others gave substantial evidence that to more accurately 
model crack turning behavior in pressurized cylinders, a 2 nd order theory was needed such 
as that described by Finnie and Saith [13]. (Here, 2 nd order refers here to the inclusion of 
the 2 nd term, or T-stress, in the asymptotic stress field in the vicinity of the crack tip, which 
is neglected in the Erdogan and Sih theory.) Knops was the first to implement this theory 
in an adaptive mesh finite element code, yet despite the improved theory his results for the 
Boeing narrow body panel test compare very closely with those of Potyondy. Like 
Potyondy, he was unable to simulate the turning of the crack in the vicinity of the stiffener, 
and the resulting flapping phenomena. Chen [10] also modeled the Boeing test panel using 
an adaptive mesh scheme and a 2 nd order theory that includes fracture orthotropy [14] 
(fracture orthotropy referring to the variation of the fracture toughness with material 
orientation, symmetric about two axes). His results for both isotropic and orthotropic 
cases are also presented in Figure 2. Note that in his isotropic run, he was for the first time 
able to simulate the turning of the crack at the tear strap. In the orthotropic case, he 
matched the first part of the trajectory considerably better, but was still unable simulate 
crack turning at the tear strap, because the crack was growing in the preferred direction of 
the 2024-T3 fuselage skin, and effect of the T-stress was insufficient to turn the crack. 
Despite this shortcoming, the ability to simulate crack deflection by a tear strap for the 
isotropic case was a significant first. 
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While there were various differences in the finite element implementation, it appears that the 
main difference in Chen’s isotropic analysis which enabled him to show turning and 
flapping where Knops did not was the characteristic length, r c , chosen for use in the 
analysis. Chen used a value of 0.09 inches which was found during finite element 
simulations to correlate fairly well with crack paths observed in double cantilever beam 
specimens reported by Pettit [4], whereas Knops used a value of 0.05 inches 1 . Attempts in 
the literature to evaluate r c often give contradictory results, and the physical meaning and 
phenomenological basis of this parameter has not been well understood (some might 
question whether such a characteristic length in fact exists at all). One object of this paper 
is to offer a physical explanation of the existence and character of the characteristic length, 
and to provide an approximate means to estimate r c from more commonly available tensile 
and fracture properties. 

Also of importance is the accuracy of the T-stress calculation. Very recently [34], a high 
polynomial order (p-version) finite element code implementation of a contour integral 
solution for computation of the T-stress was utilized to obtain T-stress values of 
unprecedented accuracy. Equally as important, an error estimation parameter was 
developed to allow a quantitative estimate of the error associated with a given crack tip 
discretization geometry (referred to hereafter as the crack tip rosette). As will be shown 
herein, much of the error is of a systematic nature, and it is possible to calibrate a given 
crack tip rosette geometry against known highly accurate solutions, thus enabling 
correction of the error in the T-stress computation a posteriori.. 


2.0 REVIEW OF 2 nd ORDER CRACK TURNING THEORY 

The mixed mode expression in polar coordinates for the elastic stress field around a crack 
tip (Figure 3) is given by (truncating after the second term of the infinite expansion given 
by Williams [15]) 
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Where Ki and Kn are the mode I and mode II stress intensity factors (Figure 4), and the T- 
stress is the constant, or far field stress component. A physical understanding of what the 

T-stress represents can be revealed by examining the radial stress of Equation (1) at 0=±r, 
for which 


1 Actually, Knops is not explicit with regard to the value of r c utilized for 2024-T3. However, he did 
quote a value for PMMA plastic plate of 1.3 mm (0.05 inches) based on the results of Ramulu and 

Kobayashi [Knps68,1 1] . Kosai, Kobayashi and Ramulu [6] later gave the same value for 2024 
aluminum, and it appears that Knops used this value for both materials. 
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In the absence of the mode II loading, the T stress is simply the stress along the crack flank 
immediately behind the crack. One could imagine that the T stress could be introduced into 
a body with a straight crack by imposing a constant far-field stress of arbitrary magnitude 
parallel to the crack without affecting the stress intensity factors. However, it is intuitively 
apparent that if this load is sufficiently increased in tension, a propagating crack may well 
want to turn and grow normal to this far field stress. 


While many crack turning theories have been proposed, one of the most promising is the 
second order maximum tangential stress theory [20,13,6], which assumes that the crack 
will propagate in the direction of maximum tangential stress occurring at a material specific 
radius, r c , from the crack tip. To maximize a e , the derivative of equation (2) is set equal to 
zero at r=r c , which can be rearranged to obtain the transcendental equation in terms of the 

instantaneous turning angle, 6 C 
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which is plotted in Figure 5 in normalized format. Note that this theory predicts turning 
(0 C < 0) in a pure mode I environment subject to the criterion 
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The common first order theory [12] can be obtained either by setting r c =0, or by neglecting 
the T-stress, which yields from Equation (5) 


K u _ -sin0 c 

(3co.fi, -1) ( 7 ) 

This method does not require evaluation of the T stress, but as apparent from Figure 4, is 
expected to lose accuracy in problems where substantial positive T stresses exist (as when a 
crack approaches a tear strap or stiffener in a pressurized fuselage). It is further noted from 
the foregoing that the effect of the T-stress vanishes unless the characteristic length is non- 
zero. 


3.0 THE CHARACTERISTIC LENGTH, r c 

3.1 Literature Review 
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Irwin [16], Dugdale [17] and others gave approximate expressions for the size of plastic 
zone in front of a crack tip in an elastic-plastic materials. Ingraffea [18] and others have 
shown that due to microcracking and strain softening, inelastic zones also appear about the 
crack tip in more brittle materials such as rock. There has been a general agreement in the 
literature since that the characteristic length associated crack turning should be smaller than 
the inelastic zones identified in these studies. 

Rice and Johnson [19] discussed the role of various characteristic lengths associated with 
microscopic failure mechanisms in elastic plastic materials, including the crack blunting 
radius and void spacing, in the context of plain strain fracture problems. 

The existence of the characteristic length associated with crack turning was proposed by 
Williams and Ewing [20], who suggested that a crack might be assumed to propagate in the 
direction of maximum tangential stress evaluated at a point a finite distance ahead of the 
crack tip. As an estimate of the characteristic length for PMMA, they referenced a previous 
work by Constable, Williams, and Culver [21] which identified equivalent flaw sizes based 
on fatigue thresholds of in polyvinyl chloride of the order of 0.0025 inches. Constable 
conjectured that the equivalent flaw effect might be associated with crazing. 

Streit and Finnie [22] determined r c using photoelastic methods for 7075-T651 aluminum 
plate to be 0.010 inches. They used the observed onset of path instability of an initially 
self-similar crack as the basis for determining r c . They described r c as the distance at which 
void growth or crack initiation will occur, referencing Rice and Johnson and others. 

Using similar methods, Ramulu and Kobayashi [23] experimentally determined r c for 
PMMA 2 to be 0.05 inches, which is quite a bit larger than had been previously expected. 
Theocaris and Andrianopolis [24] independently obtained similar results. Kosai and 
Kobayashi [25] later described r c for PMMA as “a material constant which specifies the 
characteristic crack tip region in which the off-axis micro-cracks are triggered and 
connected to the main lead crack tip”, similar to the void growth spacing assertion 
mentioned previously. 

Kosai, Kobayashi, and Ramulu [6] later estimated r c for 2024-T3 and 7075-T6 sheet to be 
0.06 inches based on the lengths of micro-crack branches observed along dynamic fracture 
surfaces of test specimens. This is considerably larger than the value given by Streit and 
Finnie for 7075-T651 plate, but the method of determination of r c is completely different 
than previous methods. Also, the reduced thickness of the sheet relative to the plate 
material would likely have played a role. 

Pettit [4], found that path instability occurred in 2024-T3 3 DCB specimens at values of r 0 at 

least as high as 0. 1 1 inches due to sensitivity of the crack path to small amounts of K„ (as 
can be observed in Figure 4 near the bifurcation at r 0 =r c ), but that the crack turning radius 
approached zero at r 0 = 0.05 inches, which was subsequently used as a (likely rather poor) 

estimate of r c . Because imperfections or perturbations giving rise to small amounts of K„ 
can be found in any real specimen [14], the onset of path instability in nominally symmetric 
specimens will always occur at an r 0 value in excess of r c , thus casting uncertainty on any r c 
values obtained in that manner. 


2 Actually, the material used by Ramulu and Kobayashi is designated as “Homelite 100” which the 
author understands is a commercial form of PMMA. Also, their experiments were dynamic in 
nature, though they claimed that static values of r c would be comparable. 

3 Actually, the material was of NASA vintage stock made to the earlier 24ST designation. 
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3.2 Present Theoretical Development 


The following development will be discussed in the context of plasticity, with the focus on 
metallic materials. Nevertheless, the general principles described can also be applied to 
materials that fracture by way of micro-cracking or other inelastic effects. 

A simple tensile test of a strain hardening material yields the familiar engineering and true 
stress-strain plots shown schematically in Figure 6 (the less familiar “Strain Softening 
Zone” shown at the end of the curve will be discussed later). In accordance with a well- 
known plastic instability theory attributed to Considere, the maximum load, F, occurs 
when the specimen rate of area reduction equals the rate of strain hardening [25] 

dF = adA + Ada = 0 (8) 


Rearranging, 
a A 


da 

— = a 
de 


(9) 


It is equally well established that the point of maximum load also defines the onset of 
localized deformation or necking in the specimen. That this true can be clearly illustrated 
by likening the specimen to a series of nonlinear springs of unit length, as illustrated in 

dF 

Figure 7. Each spring may be considered to have a local spring constant — . As the 

de 

series is stretched, all the springs elongate in proportion to their compliance, — . As the 

dF 

stiffness of any one of the springs becomes zero, then its compliance becomes infinite, all 
the other springs unload, and only that spring elongates. All along the specimen a 
compliance of zero is approached as the maximum load is approached, but due to some 
imperfection one segment reaches that point first, and necking begins there. 

Once localized deformation has begun, the location of the future failure of the specimen has 
been determined. As fracture develops, the processes which occur after the onset of 
localized deformation may differ from material to material, but the location of fracture is set 
in a macroscopic sense at the onset of plastic instability. 

Assuming strain hardening of the exponential form, 


a = ke m (10) 

and substituting this into equation (9), one obtains the true plastic strain at the engineering 
ultimate stress 



Noting that S = cxe e (where e is the base of the natural log) we obtain from (10) and (11) 
the engineering ultimate stress of the material in terms of k and m 

S u j t = krn m e~ m (12) 

Defining a o as the 0.2 percent offset yield strength, we have k=c o (,002)' m and 


q 0 J. 002c\ m 

S u i, V m J 


(13) 


Swift [26] developed similar relationships for plastic instability and necking in sheet 
material under tensile plane stress conditions. Necking is also observed in front of the 
crack tip, and as will be shown, may play a significant role in crack path formation in sheet 
metal. 


A crack growing in a thin sheet is illustrated in Figure 8, with the necking region shown in 
the vicinity of the crack tip well within the bounds of the plastic zone, since necking must 
occur after some plastic deformation as in the tensile test. It is further asserted based on 
observation that the crack will eventually develop along the necking line, and that the future 
crack path is therefore known out to the onset of necking (barring some abrupt change in 
the load environment). In order to support this notion, it is observed that the sectional load 
(load/in) distribution ahead of a propagating crack and normal to the future crack path 

(shown as a n in Figure 8) must have a maximum a finite distance away from the crack tip 
as shown schematically in Figure 9. 

McMeeking [27], using nonlinear finite element computations, has shown this to be true 
for the self-similar plane strain stationary crack due to crack blunting. The sectional load 
distribution of Figure 9 is also supported by the argument that all real materials must exhibit 
strain softening behavior across any real failure interface down to zero load 4 as asserted in 
Figure 6. 

It is now observed that in a manner equivalent to the Considere criterion mentioned 
previously, that the onset of localized deformation coincides with the instability point 
defined by 

dN n =0 (14) 


where N n is the sectional load normal to the future crack path. For elastic plastic strain 
hardening materials, this point marks the onset of plastic instability, whether in plane strain 
or plain stress. The path of localized deformation marks the future crack path just like 


4 Strain softening might alternately be described as the advanced stages of deformation 
localization, including void growth and coalescence for metallic materials. The interface has finite 
residual strength until the last two atoms of a given interface separate, (and even as they separate, 
they do so with a smooth load/displacement relationship). Given that the atomic bonds, even in a 
tensile test, must break in some sequence, and cannot separate at exactly the same moment, it 
follows that in the limit of absolute displacement control the failure of any interface could be 
defined as a quasi-static progression of damage as the load drops smoothly (at some scale) to 
zero. The failure interface of a slow stable tearing interface approaches this limit of absolute 
displacement control [18] 
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necking in the tensile specimen predetermines the eventual failure location. The instability 
point may thus be considered as the end of the known future crack path, the point where 
the material is “deciding” where the crack will go next. Presumably as the crack grows the 

instability point would migrate to the location where N n is maximized. Thus the distance 
from the physical crack tip to the instability point, l c , might be considered to physically 
represent the characteristic length associated with crack turning. 


Unfortunately, evaluation of the generally curvilinear l c is a daunting task. For the 
purposes of this study, it shall be assumed to be straight, and of approximately constant 
length, r c for a given material and thickness. In order to obtain an estimate of that length, 
the simplest case of a self-similar crack will suffice. For the plane strain case, McMeeking 
and Parks [28] have already shown that the maximum stress occurs for materials with 
moderate strain hardening at a distance approximately given by 



(15) 


Where for small scale yielding and mode I loading, the strain energy release rate J, is 
related to the stress intensity factor by 


J 



(16) 


Thus for plain strain we may write 


jsL 

<7 0 E 


(17) 


It should be cautioned, however, that while this expression is a reasonable estimate for 
plain strain with T=0, it has been shown [40] that the distance to the maximum stress point 
varies with the T-stress, particularly if it is negative. 

For the plane stress case, thickness strain (necking) is possible. We can thus rewrite (14) 
as 


dN n - dNfl = tdcr q + CTfldt = 0 

(18) 


(19) 

Oq t 


Assuming incompressibility. 


-de z = de 0 (l + p) 

(20) 


where p = 


h. 

£e 


Using the Von Mises yield condition, the differential equivalent strain can be expressed as 
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( 21 ) 


dc — d £q ^ {\ + p + p j 


Combining (19), (20), and (21) and noting that for a constant p, the equivalent stress 

da d era 

conforms to the equality = we obtain 

a a e 


da aV3(l + p) 
de 2(l + p + p 2 ) 


( 22 ) 


Assuming an exponential strain hardening relationship for the equivalent stress and strain 
similar to (10), we obtain the critical strain and stress at the instability point. 


2m|l + p + p 2 j 
;crit= V3(l + p) 


(23) 


^crit 


= k 


2m(l + p + p j 
V3(l + p) 


-im 


(24) 


Assuming the necking region in plane stress is somewhat larger than the crack blunting 
effected zone in plane strain (which will be observed hereafter), we can obtain an estimate 

of p based on the Hutchinson Rice Rosengren (HRR) asymptotic field for plane stress 
cracks in strain hardening materials [29,30, 31]. This inherently involves the assumption 
of proportional flow, which for a propagating crack is less than desirable. Directly ahead 
of a straight propagating crack, the flow should be fairly proportional, but the deviation 
from proportionality will have an unknown effect on the results forward of the crack. Also 
of concern is the fact that the HRR field is derived based on the assumption of small strain 
theory. The very existence of plastic instability is by definition a large strain effect. 
Nevertheless, McMeeking and Parks did observe in their analysis that the plane strain HRR 
field was valid up to the point of maximum stress. Note also that true stress and strain and 
engineering stress and strain are fairly close up until plastic instability for most structural 
materials. Proceeding with the above cautions in mind, the effective stress given by the 
HRR solution, written here in terms of the far field stress intensity factor, is 

1 

n + 1 

a e (n,0) (25) 


a = a r 


K? 


G^o I n r . 


For 0 = 0, Hutchinson normalized the a e term to unity for the plane stress case. The 
Ramberg-Osgood material parameters a and n are related to the exponential strain 
hardening parameters from Equation (10) by 


k = a ( / 1 1/n ) 



(26) 
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m = 1/n 


(27) 


Hutchinson gives numerical results for I n and p which are functions of the exponent n. 

For the propagating crack, we may assume that K, is equal to the propagating value, which 
for stable tearing may be denoted K c . We may thus approximately equate (24) and (25) at 
r=r c , and to obtain after some rearrangement 


r c ~ 


^V3 cr 0 an(l + p) ' 

(n+1) 

n 

’ K 2 

2 E(l + p + p 2 )J 


aa 2 ! n _ 
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Equation (28) can be combined with equations (12), (26), and (27) to obtain an expression 
in terms of the engineering ultimate strength 


(i + p) > 
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For n<7, p makes little contribution, and I n is nearly linear at least up to the maximum 
value of n=13 given by Hutchinson, and probably well beyond. Thus, in this range we 
can further approximate with no significant loss of accuracy 
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(30) 


If a suitable value for n is not available, one can obtain an approximation using the ratio of 
yield and ultimate strength from the implicit equation 


-^- = (.002ne) 1/n (31) 

Suit 

which was obtained by combining equations (13) and (27). 

A predicted r c value for 2024-T3 alloy is given in Table 1 based on Equations (30) and 
(31). The mechanical properties are B basis values from MIL-HDBK-5G [32], and the 
fracture toughness the maximum R-curve value obtained based on 48 inch wide R-curve 

panels reported by Gruber et al [33]. The predicted r c value is in satisfactory agreement 
with the value chosen by Chen to match observed crack paths. One caution, however, is 
that the value of Kc used is only reached after 4-5 inches of stable tearing, whereas the 
coupon test data achieved crack turning after less than an inch of growth, at which point the 
characteristic length associated with necking should have been considerably smaller. 


Table 1. Calculation of r r for 2024-T3 


.063 inch, 2024-T3 clad 

r c 

■ 

H 

S u .t 

(ksi) 

n 

K c 

(ksiVin) 

E 

(ksi) 

K c 2 /S UIt E 

Plane Stress 

LT 

42 

62 

8.04 

180 

10500 

0.050 

0.098 
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4.0 CALCULATION OF ACCURATE T-STRESS VALUES 


4. 1 Literature Review 


T-stress calculations have been performed by various authors. Larsson and Carlsson [35] 
evaluated the T-stress using finite elements. Leevers and Radon [36] directly imposed the 
infinite series solution given by Williams [15] in a variational approach to obtain estimates 
of K, and T simultaneously. They gave estimates of the T-stress in the form of the 
dimensionless parameter B, where 


B = 


TVa 

Ki 


(32) 


Based on the convergence observed, Leevers and Radon estimated the error in the B values 
provided to be less than three percent. Sham [37] used 2 nd order weight functions and a 
work conjugate integral to calculate T-stresses in various specimen configurations. Fett 
[38, 39] introduced a Green’s function approach to calculate T-stresses, and analyzed 
numerous configurations. A more approximate displacement correlation method was 
outlined by Al-Ani and Hancock [40] which is nevertheless easy to implement in plate and 
shell codes, and has been utilized in various forms by other authors [4, 9, 10]. 

Cardew et al [41] and Kfouri [42] computed the T-stress using a modified J-integral based 
on unpublished work of Eshelby, and also gave results for selected specimens based on 
finite element analyses. Another type of path independent integral based on the Betti- 
Rayleigh reciprocal theorem has also been proposed [43,44], and shown to be 
mathematically equivalent to the J-integral method by Chen et al [34]. By implementing the 
contour integral solution into a high polynomial order (p-version) finite element program, 
Chen obtained T-stresses that were claimed to be numerically exact to six significant 
figures. The numerical accuracy was verified by way of an exact benchmark solution (a 
crack tip and surrounding region with the exact boundary conditions applied corresponding 
to arbitrary combinations of K,, K„, and T) and a theoretical error relationship 


HP r T~' — ^ I 

e T = t FE -T = e T -J= 

Jr, 


(33) 


where e T is the error in the computed T-stress, r, is the characteristic dimension of the 
integration zone (typically taken as external radius for circular zones, half side length for 
square zones), and e T represents the discretization error in the vicinity of the integration 
zone. 


4.2 Error Correction 


Equation (33) was derived by recognizing that the stress contribution of the singular terms 
in the stress field will converge far slower than the contribution of the non-singular terms, 
leaving an error in the coefficients of all terms proportional to the coefficients of the 
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singular terms 5 . The square root term in the denominator was included due to dimensional 
considerations, consistent with the form of Equations (1-3). Based on the convergence rate 
argument, terms of higher order than T are expected to contribute little error, though this 
was not checked in the benchmark example, since higher order terms were excluded. 

Note that unlike the error estimation expression given in [34], Equation (33) is given with 
no absolute value signs to enforce that the error measure always be positive. As shown in 
Figure 10 for the integration path and rosette geometry of Chen, ej is virtually constant 
and characteristically negative 6 for a element polynomial order, p, above 12. The error is 
thus systematic, and can thus be corrected as shown in the next section. The relative error 
of the T-stress can be expressed in a manner consistent with the characteristic sign of the 
error. 


e T _ - Kj 

eT '« = m T mfi 


(34) 


The assertion that the relative error scales with K/Wr, is supported by the fact that 
geometrically similar finite element models that differ only by a scaling factor (which also 
implies that the integration zone is likewise scaled) will have the same relative error. 
Considering the rosette as a finite element model with imposed boundary conditions 
representing K, and T, and recognizing the similarity of all K, and T fields relative to a 
characteristic length (K,/T) 2 , one may therefore conclude that the combination of such a 
field with a rosette model of fixed geometry and scale relative to the field characteristic 
length will be similar (and thus have comparable relative error) to all other rosette/field 
combinations with the same relative scaling ratio. 


Because e T represents the discretization error in the vicinity of the integration zone, it 
should thus be relatively constant so long as the mesh geometry, or rosette, within the 
integration zone is geometrically similar for all problems. The mesh geometry outside of 
the integration zone is of secondary influence, and may change from problem to problem, 
thus its effect will be treated as a probabilistic source of error. Nevertheless, provided that 
the external mesh is reasonably proportioned, the error introduced should be relatively 
small. 


Not noted in the previous paper is the fact that the err or in K, likewise scales to itself, 
e Ki =k Ife -K I = ^/ K I 

No characteristic length is necessary because all K fields are similar without respect to a 
characteristic length. Thus, the relative error in K, is approximately constant for a given 
integration path/rosette geometry 


5 K„ was also found to pollute T slightly, but at a level orders of magnitude lower than K,. It would 
thus likely have even a lesser effect on K,, thus it is hereafter neglected. In nearly pure mode II 
loading, however, the effect of K,, might be significant, and could be normalized in a similar 
manner. 

6 Note that the characteristic sign consistency of the error may not hold if the crack tip elements 
are included in the integration zone, as will be explained. Also the sign of the error may change 
with p. 
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( 36 ) 


e K, - 

eK Irel = ~^~ €K 1 

As shown in Figure 1 1, e^is constant and characteristically positive for the integration 
path/crack tip rosette geometry of Chen p>9, after which the benefit of increasing p is lost 
for this rosette geometry, and error of the opposite sign appears to take over. A similar 
correction is surely possible for K„, but is not discussed at present. 

If the sign of the error is always the same. Equations (33) and (36) can be utilized to correct 
the finite element solutions for K, and T. Designating T adj and K ladj as the adjusted values 

T ~ Tadj = T FE ' e T ~j^ 


K I ~ ^Iadj “ ^IFE ~*K[K I 

= k ife/( 1+ ^/) 


(38) 


The adjusted values so obtained will still vary from the exact, because in practice e T and 
e Kl are distribution functions, and can be characterized with mean values and standard 

deviations. As mentioned previously, the random component of the error may be 
associated with mesh variations outside the rosette and solution oscillation at the crack tip. 
In fact, if the equivalent domain integrals for K, and T are evaluated including in the domain 
the crack tip elements themselves, the random component of error may be larger than the 
systematic error, resulting in overall errors which may vary in sign, and are thus less 
correctable. For this reason it is recommended that the crack tip elements be excluded from 
the integration zone. 

In order to evaluate Equations (37) and (38), one must have an estimate of the calibration 
factors e T and e Kj . For individual cases with known exact or highly accurate reference 

solutions, sample values of ej and e^ / can be obtained as 


= (Tfe ■ T)/ (39) 

^K, = e Kj =TfE - T ( 40 ) 

The mean and standard deviation can estimated based on models of a sampling of reference 
solutions, preferably over a range of stress intensity factors and integration radii. For the 
data in Figures 10 and 11, three samples of each error parameter were evaluated at each 
polynomial order p. Average values of the correction factors are given in Table 2 given to 
eight computed places, which are likely all significant for p=l, but perhaps only one digit is 
significant at p=12. 
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Table 2. Average Values of the Correction Factors 





p 

e K, | 

e T 

1 

4.7282230E-01 

-1.8830976E-01 

2 

1. 006327 IE-01 

-2.73707 18E-02 

3 

2.1770302E-02 ~1 

-5.8434480E-03 

4 

4.5650887E-03 

-6.7037228E-04 

5 

9.1459033E-04 

-1.9170107E-04 

6 

1.7818900E-04 

-5.0844804E-05 

7 

3.2591667E-05 

-7.4061435E-06 

8 

4.3786667E-06 

-2.3991630E-06 

9 

-3.1936667E-07 

-4.2137896E-07 

10 

-9.3926667E-07 

-9.3609849E-08 

11 

-8.1720000E-07 

-2.1412723E-08 

12 

-6.4710000E-07 

2.5544523E-09 


Using these values, the solutions were corrected resulting in the miniscule errors in T adj and 
KI ad indicated by the error bands in Figures 10 and 1 1. Note that the error obtained in this 
manner is independent of p. Nevertheless, this example is not realistic because the models 
were completely similar — not just similar in the vicinity of the integration zone. 

In order to obtain a more realistic representation, K, and T estimates from the DCB 
specimen geometry presented in [34] was adjusted according to Equations (37) and (38) to 
obtain the data given in Figure 12. Each analysis used the same six-tiered rosette used in 
the benchmark example, with the integration zone encompassing the first layer of elements 
outside the crack tip elements as before, thus the correction factors given in Table 2 were 
applicable. As expected, the benefit of correction in this real case is nearly as good as for 
the benchmark solution. Nevertheless, the improvement in the accuracy of KI is typically 
an order of magnitude or more, though not as good as at low p values. The improvement 
in T due to correction is about half an order of magnitude. 

The rosette geometry utilized, however, is very poor at low p values, which are common in 
available finite element codes. FRANC2D, a two-dimensional fracture code developed at 
Cornell University, utilizes a more suitable crack tip rosette made up of eight six-node 
quarter-point singular elements as shown in Figure 14. FRANC2D is being used as a 
platform to develop crack turning analyses, thus it is of interest to calibrate the rosette for 
T-stress calculations. Due to the use of singular elements, the accuracy of stress intensity 
factors calculated by FRANC2D is typically within one percent, thus no correction is 
necessary there. 

Figure 14 shows a rosette subdivided once with the inner radius equal to one-half the outer 
radius. The standard version of FRANC2D allows the user to manipulate the number of 
subdivisions and the subdivision ratio. A subdivision ratio of 0.5 was selected for this 
study. Also, the version of FRANC2D utilized has been updated to include the contour 
integral solution for T. It was also found that FRANC2D models cracks with a small initial 
crack opening, which causes numerical difficulties when evaluating the contour integrals. 
As a temporary solution, gap elements with zero interference were specified to close the 
crack prior to running the analyses. 
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In order to provide data from which to evaluated , several runs were made, analyzing 
selected specimens for which highly accurate reference solutions were obtained by Chen et 
al. The error in the computed T-stresses is shown in normalized fashion in Figure 15. For 
this rosette geometry/integration path of Figure 14, it was found that 

ej = -0.00825 ± .00255 

(Mean) (Std. Deviation) 

A plot of the corrected data with lines denoting various confidence levels is given in Figure 
15. The average error (50% confidence level) of the corrected solution was about on fifth 
of the original error. 


5.0 CONCLUSION 

The characteristic length associated with crack turning, r c , has been characterized as the 
distance ahead of the crack tip at which plastic instability occurs for both plane strain and 

plane stress fracture problems in strain hardening materials. An estimated plane stress r c 
value has been determined for 2024-T3, and compares favorably with a recent value given 
in the literature which results in good correlation with test data. Nevertheless, due to the 
many severe assumptions associated with the plane stress solution caution is recommended 
with regard to its use pending further evaluation. 

An error correction methodology has been developed for the T-stress, which can also be 
applied to stress intensity factors. With this methodology, the accuracy of a contour 
integral based T-stress algorithm has been substantially increased, and an estimate of the 
remaining error is provided. 
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Figure 1. Crack Turning and Flapping in Boeing 707 Test Panel [1] 
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— o— predicted (Chen, revised boundary condition) 

— a— predicted (Chen, Potyondy boundary condition) 
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a.) first order maximum tangential stress theory 


— □ — measured 

-o— predicted (Chen, Isotropic, rc=0.09 in) 

— predicted (Chen, Oithotropic, rc=0.09 in) 



b) Second order maximum tangential stress theory 


Figure 2. Correlation of Crack Growth Simulations of Various Investigators (from [10]) 
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Figure 3. Polar Coordinate System About a Crack Tip, also Showing Physical Meaning of 

the T-Stress 



Figure 4. Mode I and Mode II Loading of a Cracked Body (Shown with Positive Sign) 
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Figure 7. Illustration of Instability Principle Using Nonlinear Springs 



B-B 


Figure 8. Schematic Physical Process Zones in the Vicinity of a Propagating Crack 
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Figure 10. Accuracy Assessment of T Computations Using p-version Contour Integral 

Implementation of Chen et al [34] 



Figure 11. Accuracy Assessment of K, Computations Using p-version J-Integral 

Implementation of Chen et al[34] 
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Figure 12. Comparison of Error in as-Calculated and Corrected Values of K, and Tina 
Double Cantilever Beam Specimen, h/w=.2, a/w=. 5 



Figure 13. Geometry and Integration Path Used for FRANC2D Rosette Calibration 
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Figure 14. As-Calculated Error in T-stress Using FRANC2D 



Figure 14. Corrected Error in T-stress Calculated by FRANC2D with Upper Bounds 

Representing Estimated Error 
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Abstract 

Elastic-plastic Crack Tip Opening Displacement (CTOD) criteria proposed for crack 
propagation and turning [1,2,3] are examined in the context of linear elastic fracture 
mechanics. In particular, the use of the CTOD as a crack turning criterion is discussed as 
an extension of the mode I dominant criterion of Cotterel and Rice [4] by simulating 
plastic deformation via a virtual crack kink which opens in pure mode I or pure mode n. 
The mode I dominant case is found to result in a crack turning theory similar to the 
second order maximum tangential stress theory with a transition to mode II deformation 
anticipated under certain conditions. 


INTRODUCTION 

Crack turning has been studied by many investigators [5], and is well recognized as a 
potentially important phenomenon in structural applications, such as crack turning and 
flapping in aircraft fuselage structures [6, 7]. Crack path prediction using linear elastic 
maximum tangential stress theory [8, 9] has been most commonly used in structural 
analysis. For crack turning and flapping phenomena and other problems, the effect of the 
second order, or T-stress term in the asymptotic stress field has been shown to be very 
important [10, 11]. Approximate theories have accounted for the directional variation of 
fracture properties in anisotropic materials [12, 13]. Inclusion of the T-stress term in the 
linear elastic theory requires the use of a material dependent characteristic length parameter, 
r c , associated with the process zone size. Because the process zone presumably sizes with 
the plastic zone, it follows that the effect of the T stress will increase with the extent of the 
plastic zone. This observation has been used to explain quantitatively why slow fatigue 
crack growth turns more gradually than static tearing of thin aluminum sheets [11]. 

As computational hardware and software have become increasingly powerful, nonlinear 
finite element analysis of complex structures has become more viable, leading to an 
increased emphasis on elastic-plastic methods to predict crack propagation and trajectory in 
situations where the effects of plasticity are potentially significant. To date these efforts 
have focussed on isotropic materials. Early work in this field showed that plasticity can 
effect the crack trajectory [14], and supported a notion that crack propagation and trajectory 
could be correlated with critical CTOD values [1, 2]. Sutton et al [3] recently proposed 
curves based on results from a small scale yielding (SSY) boundary layer crack tip model 
and laboratory experiments which can be used to infer the crack growth direction from the 
ratio of the mode I and mode II CTOD components. The curves include a transition from 
mode I to mode II dominated crack extension that has been observed by previous 
investigators [2, 15], and were used to predict stable tearing crack paths in various 
specimens, with very encouraging results. 
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The CTOD crack trajectory curves appeared to be independent of the T-stress [3], and yet 
seemed to adequately capture the small turning radius of statically tearing cracks in a high T 
stress environment. Thus, questions were raised about the influence of the T-stress on 
crack trajectory, and how or if its effect is captured by the CTOD approach. In an attempt 
to address these issues, the CTOD method is examined from a surrogate linear elastic 
standpoint, once again including the effects of plasticity and the T-stress by virtue of a 
characteristic length. 


2 nd order linear elastic surrogate for elastic-plastic theory 

Consider a lead crack under plane stress conditions with a plastic zone as shown in Figure 
la. Compared to an elastic crack, the plastic zone results in additional deformation that can 
be approximated by a virtual elastic crack kink as shown in Figure lb. For self-similar 
crack growth, Wells [15], used the Irwin plastic zone correction as an approximation of the 
effective elastic kink length to obtain an estimate of the CTOD. While the appropriate 
choice of length may remain in question, it is not unreasonable to assert that for a given 

material and loading, there is a unique kink length, b c , and orientation, 0 C , which will best 
simulate the deformation field as one moves away from the crack tip into the elastic 
region 1 . It is further postulated that a crack kink so defined would provide a reasonable 
approximation of the future crack trajectory. For a crack propagating under steady-state 
conditions, b c would be expected to assume a constant, material-dependent value. 


As the case of perfect elasticity is approached (as for so-called brittle materials, and also 
approximately for slow fatigue crack growth), the length of the virtual kink necessarily 
vanishes. In this limiting case, Cotterel and Rice [4] concluded that the crack propagates in 
pure mode I, which is equivalent to the criterion K„= 0. For the finite virtual kink, it is not 
clear that this will continue to be the case, and based on the experimental observations [2, 
15], depending on the loading conditions, cracks are observed to develop trajectories 
corresponding either to pure mode I or pure mode II crack opening. For the time being, it 
will be assumed that the virtual kink tip will operate at pure mode I, but the possibility of 
pure mode 13 propagation will also be considered. 

Isida and Nishino [17] give a solutions for a crack in an infinite plate with a kink at one end 
which can be superimposed to obtain results for general K,, K n , and T loading. The stress 
intensity factors at the kink tip, k,, and k n , may be expressed as 

k„ = f$>K, + F}PK„ - FjPTjH 

Where a is the crack length, and Fj° are functions of the kink angle, 9, and the normalized 
kink length, b/a. The crack length parameter can be eliminated by normalizing in the form 


1 It is suggested for future study that the virtual crack kink dimensions be selected to give an 
equivalent J integral (evaluated at a radius somewhat large compared to the kink length) to the 
elastic-plastic crack. 
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where -Jb/a factors cleanly out of functions F„ (2) , and T is defined with b=b c in a manner 
similar to that utilized for the maximum stress theory [1 lpettit97]. 


T = -—^2nb c 
3 K/ c 


(3) 


Values of the crack propagation angle, 6 C , can be obtained by varying 6 to enforce K„ = 0 

for various combinations of K,/K„ and F, as presented in Figure 2a. Also shown for 
comparison is the 2 nd order maximum tangential stress theory in normalized format [9, 7 , 
11] with its characteristic length, r c . Note that the two theories are nearly equivalent 
(though not identical) if one recognizes that the characteristic lengths differ at the 
bifurcation by a constant factor, 

b c = 2.2 lr c (4) 

Considering that the 2 nd order maximumum tangential stress theory is well regarded for its 
ability to predict physically observed phenoma, this similarity is encouraging. For the ideal 
linear elastic case, b= 0, and the crack turning criterion becomes a function of K,/K, 
independent of the T-stress. This is plotted in Figure 3 in terms of the crack tip opening 
displacements using the notation of Sutton [3] 


a = arctan 


Bil 

D,) 


(5) 


Where D , and D „ are the mode I and mode II components of crack opening, with the 
observation that for the linear elastic case, D,/D, = K,/K r 

For the kinked-crack representation of the elastic-plastic crack tip, the crack tip opening 
displacements may be approximated as indicated in Figures 4a and 4b, depending on 
whether the virtual crack is assumed to develop in pure mode I or pure mode II. Assuming 
the CTOD is measured in a finite element simulation some small distance, d, behind the 
physical crack tip (the physical crack being represented by the base of the kink), the 
following relations may be obtained. 


-sin0 c 

for a pure mode I kink 

(6a) 

cos 6 c +d/b c 

~ - tan 0 C 

(approximation for d/b c -$ 0) 

(6b) 

II 

O 

o 

1—*- 

for a pure mode II kink 

(7) 
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Clearly these equations are approximate, but are most accurate and meaningful for small 
d/b c . As d/b c becomes large, D t / D , must approach the elastic K,/K, of the physical crack 
tip, and the effect of plasticity and the T-stress on the turning angle is lost. In the limit as 
d/b c ->0, the Equations (6,7) simply state that for mode I dominant growth, the crack grows 
perpendicular to the CTOD, and for mode II dominant growth, the crack grows parallel to 
the CTOD. Based on the linear elastic infinitesimal kink theory, no combination K, and K„ 
at the tip of the unkinked crack results in a pure mode I kink at an angle above 74.45 
degrees. The theory predicts this limiting value only when the lead crack is loaded in pure 
mode II (the corresponding value predicted by the maximum tangential stress theory is 70.5 
degrees). With a finite (virtual) kink length, the limiting value becomes a slight function of 
the T-stress, but only at very long kinks or very high T-stresses. If, as postulated for 
isotropic materials, the virtual kink must be pure mode I or pure mode II, then one might 
expect that the crack would transiton to a non-mode I failure at kink angles approaching 
74.45 degrees. Sutton observed a trasition to mode II failure at kink angles of about 70 
degrees in 2024-T3 aluminum, under partially mode I loading conditions. Note however, 
there appears to be no transition whatsoever in materials such as plexiglass, which under 
pure mode II loading have been observed to kink consistently at angles very near the 
maximum stress theory value of 70.5 degrees [18]. 

Sutton [3] reported 2-D simulations 2 using an elastic/small-deformation plasticity model of 
a edge crack in a 30.4 inch radius circular plate of 2024-T3 aluminum. The crack tip was 
centered on the model, and traction boundary conditions were applied representing various 
amounts of K,, K„ and T. For a given loading, the CTOD components D, and D u were first 
measured at a distance <7=0.04 inches behind the crack tip, after which the crack was 
extended 0.0076 inches as a physical kink in several trial directions to find the kink angle, 

0 C that resulted in the maximum total CTOD a small distance back from the kink tip. The 
total CTOD in this case is given by 

CTOD tM = Jsf+ti, (8) 

where 8, and 8 n are the CTOD components near the kink tip. Sutton found that the kink 
orientation that produced the maximum CTOD also produced nearly pure mode I or mode II 
displacement near the kink tip, which further supports the foregoing assumptions with 
regard to the mode mixity of the virtual kink. 

With the plane stress plastic zone for 2024 on the order of inches, it is expected that d/b c for 
Sutton’s model is sufficiently small to allow the tangent approximation of Equation (6b). 
Equations (6b) and (7) are compared with test data and CTOD-based crack turning angle 
data Figure 3, showing good correllation for both mode I and mode II cases. 


DISCUSSION 

Note that the development of Equations (6) and (7) makes no assumption with regard to the 
T-stress, or even K, or K ir The primary assumption is that the opening of the virtual kink 
will be either pure mode I or pure mode II. Inasmuch as Figure 2a is also based on pure 


2 Several of the details given here were not included in the referenced paper, but were provided 
via private correspondence with Dr. Sutton. 
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mode I opening of the virtual kink tip, it is clear that several different combinations of K,, 
K„, and T acting about the physical crack tip can result in the same kink angle. Thus, it 
appears that D,/D, measured near an elastic-plastic crack tip will be fairly uniquely related 

to 0 C , but will not be uniquely related to K,/K, if plasticity is significant, and the T-stress is 
substantially non-zero. 

Also worth mentioning is the likelihood that even as the crack tip nears a stiffener or other 
geometric feature, the virtual kink might still open in pure mode I or mode II. While the 
plastic zone shape would be distorted by the presence of such irregularities (which would 
likely affect the crack turning angle), the assumptions associated with the derivation of 
equations (6) and (7) are still valid with sufficiently small d, and the crack trajectory should 
remain a substantially unique function of D,/ D , . 

It is significant to observe that the driving features of the crack turning phenomenon appear 
to be adequately represented by a small deformation theory finite-element method used by 
Sutton. In particular, such a method neglects near-tip strain localization such as necking or 
crack tip blunting, and void growth and coalescence. It thus appears that the length scale 
associated with crack trajectory is larger than the length scales associated with these near-tip 
phenomena. 
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a) Physical crack tip and plastic zone 



Virtual 

Kink 


b) Physical crack with virtual kink 

Figure 1. Crack Tip Plastic Zone Deformation Simulated by an Effective Virtual Kink 


8 








Figure 3. Comparison of Various CTOD related Crack-Turning Curves, including the data 

of Sutton [3] 
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CTOD 

measurement 

location 


a). Pure Mode I 




b) Pure Mode II 

Figure 4. Approximate Method of Obtaining CTOD Mode Mixity for Pure Mode I and 

Mode II Virtual Crack Kinks 


11 



